Fast method for dynamic MR imaging

ABSTRACT

Generalized series-based image reconstruction as used in dynamic imaging for high-speed imaging with limited k-space coverage for each time frame. Further, in acquiring low resolution data for a plurality of image frames, a full k-space data set is generated for each time frame with the measured low-resolution data and high spatial frequency data generated by the GS model constructed based on the high-resolution image(s). The algorithms of the invention have computational complexity of O(N log N) and are capable of producing high-resolution dynamic images with a small number of Fourier transform samples.

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY SPONSOREDRESEARCH OR DEVELOPMENT

[0001] The U.S. Government has rights in the disclosed inventionpursuant to NIH Grant No. 1-R21-HL62336-01 and NIH Grant No.P41-RR09784.

CROSS-REFERENCES TO RELATED APPLICATIONS

[0002] NOT APPLICABLE.

REFERENCE TO A “SEQUENCE LISTING,” A TABLE, OR A COMPUTER PROGRAMLISTING APPENDIX SUBMITTED ON A COMPACT DISK

[0003] NOT APPLICABLE.

BACKGROUND OF THE INVENTION

[0004] This invention relates generally to dynamic imaging, for example,using magnetic resonance (MR) techniques, and more particularly theinvention relates to a fast method of image reconstruction using ageneralized series for recovering unmeasured data in imagereconstruction to enable high-speed imaging. The invention findsparticular applications in imaging of time-varying object functions, forexample, following wash-in/wash-out process of injected contrast agents,monitoring perfusion/diffusion changes in brain functional studies, andmeasuring relaxation times in multiple T1-weighted or T2-weightedexperiments. However, the invention is also applicable to other imagingmodalities including X-ray CT (computerized tomography), PET (positronemission tomography), and the like.

[0005] Many MR imaging applications require collecting a time series ofimages. Conventionally, these images are acquired independently, and asa result often compromise spatial resolution for high temporalresolution or short imaging time. To overcome this problem, severaldata-sharing methods have been proposed to improve imaging efficiencyand speed. One such method relies on data sharing from one temporalframe to another, such as reduced encoded imaging by generalized seriesreconstruction or RIGR. A unique feature of this method is that one ormore high resolution reference images and a sequence of reduced dynamicdata sets (usually in central k-space) are collected. Assuming that NFourier encodings are collected for the reference data set(s) and Mencodings are collected for each of the dynamic data sets, a factor N/Mimprovement in temporal resolution is gained with this data acquisitionscheme as compared to the convention full scan imaging method.

[0006] One known method, referred to as keyhole, fills in unmeasuredhigh frequency data simply with the data from the data reference set orsets. The RIGR method recovers unmeasured data using the generalizedseries (GS) model, of which the basis functions are constructed based onthe reference image(s). A key step in GS model-base image reconstructionis the determination of the series coefficients, which involve solving aset of linear equations. Although the computational problem ismanageable for k-space data truncated only along one direction, it canpresent a significant problem for data sets truncated in multipledirections as often encountered in multidimensional imaging. To addressthis problem, a fast approximate solution has been proposed forapplication in 3D time-resolve angiography. See Madore and Pelc, “A newway to perform 3D time-resolve angiography,” Proceedings of the EighthAnnual Meeting of ISMRM, Denver, 2000; p. 697.

[0007] The present invention provides a fast algorithm forreconstructing GS model-based images. The algorithm is stable, dataconsistent, and capable of capturing large signal variations withrespect to the reference image(s):

BRIEF SUMMARY OF THE INVENTION

[0008] In accordance with the RIGR imaging scheme, one or morehigh-resolution reference images and a sequence of reduced dynamic datasets are collected.

[0009] In accordance with the basic generalized series (GS) model, adynamic image I(x) is expressed as a product of a weighting function(determined from one or more high resolution reference images) and alocal contrast modulation function, which is updated from one imageframe to another. The contrast modulation function includes seriescoefficients which are determined by solving a set of linear equations.

[0010] In accordance with the invention, a fast algorithm finds anapproximate solution for the sets of linear equations which are thenused in GS model-based image reconstruction.

[0011] In accordance with another feature of the invention, theweighting function of the GS model is determined from the magnitude ofthe reference image(s), and a regularization constant is added to theweighting function to provide stability but without adversely affectingaccuracy.

[0012] In accordance with yet another feature of the invention, thecontrast modulation function is determined from the filtered Fouriertransform of the dynamic data sets and synthetic Fourier data of theweighting function.

[0013] In accordance with yet another feature of the invention, thephase information of the reference images and measured dynamic data areutilized in every time frame along with the estimated data to therebyprovide consistency of the reconstructed image.

[0014] The invention and objects and features thereof will be morereadily apparent from the following detailed description and appendedclaims when taken with the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

[0015]FIG. 1 shows the k-space sampling pattern of the RIGR imagingtechnique involving one or more high-reference data sets. The presentinvention can also handle k-space data sets collected in other patterns,such as radial, and spiral trajectories.

[0016] FIGS. 2(a)-1(e) are images attained in a contrast-enhancedangiographic image experiment including (a) pre-injection andpost-injection reference images; (b) high resolution post-injectionFourier images reconstructed from 144 phase encodings; (c) keyholereconstructions from 16 phase encodings using the second image in (a) asthe reference; (d) GS model-based reconstructions from the same data setin (c) using the fast algorithm of the present invention; and (e) GSmodel-based reconstruction using prior art TRIGR (two referencereduced-encoding imaging by generalized-series reconstruction)algorithm.

[0017]FIG. 3 is a graph of reconstruction errors normalized against theenergy of high-resolution Fourier images for Keyhole and for Fast andOriginal TRIGR.

[0018] FIGS. 4(a)-3(e) are diffusion-weighted images acquired with avariable diffusion-weighting gradient: (a) Two high-resolution referenceimages; (b) High-resolution Fourier images reconstructed from 128 phaseencodings (horizontal direction); (c) Keyhole reconstructions from 16phase encodings using the second image in (a) as the reference; (d) GSmodel-based reconstructions from the same data set in (c) using the fastTRIGR+ algorithm of the invention, and (e) GS model-basedreconstructions using the original TRIGR+ algorithm.

[0019]FIG. 5 is a graph illustrating reconstruction errors normalizedagainst the energy of high-resolution Fourier images for Keyhole, FastTRIGR, and Original TRIGR.

DETAILED DESCRIPTION OF THE INVENTION

[0020] The basic aspects of magnetic resonance imaging (MRI), or imagingusing MRI signals, are by now well understood. Spatial information isencoded into the acquired MRI signals using magnetic field gradientsthat are controlled by a pulse sequence program. Thus, each acquiredsignal is associated with a spatial encoding parameter. Generally, eachacquired data point represents a sample of the Fourier transform of theimaged object. The function “space” in which the Fourier transformsamples are defined is called “k-space.” Often, the k-space samples areon a Cartesian grid and each pulse sequence repetition collects sampleson one line in k-space. An important spatial encoding parameter in thiscase is the so-called phase encoding amplitude. In another MRI method,the samples collected in one pulse sequence repetition fall on spiralpaths in k-space, and k-space is filled using multiple rotations of thisspiral path. An important spatial encoding parameter in this case is theangular rotation of the spiral path.

[0021] Many imaging applications involve acquiring a time series ofimages. To improve imaging speed, several “data sharing” methods havebeen proposed which collect one or a few high resolution references anda sequence of reduced data sets. In image reconstruction, two methods,known as “keyhole” and RIGR (reduced encoding imaging by generalizedseries reconstruction), have been used. Keyhole fills in the unmeasuredhigh frequency data simply with those from the reference data sets,whereas RIGR recovers the unmeasured data using a generalized seriesmodel, of which the basic functions are constructed based on thereference image or images. The present invention provides a fastalgorithm for generalized series-based image reconstruction.

[0022] The basic GS model can be written in the following form:$\begin{matrix}{{I\left( \overset{\rightarrow}{x} \right)} = {{w\left( \overset{\rightarrow}{x} \right)}\quad {\sum\limits_{n}{c_{n}^{{- }\quad 2\quad \pi \quad {{\overset{\rightarrow}{k}}_{n} \cdot \overset{\rightarrow}{x}}}}}}} & (1)\end{matrix}$

[0023] where the weighting function w({right arrow over (x)}) is chosento absorb any available a priori information, and the seriescoefficients c_(n) are chosen to assure that I({right arrow over (x)})is data-consistent. For example, w({right arrow over (x)}) could be ahigh-resolution reference image and the second term provides a localcontrast modulation. The optimality of the model can be justified fromthe principle of minimum cross-entropy, and some of its desirableproperties have been discussed in the context of constrained imagereconstruction.

[0024] A key step in GS model-based image reconstruction is thedetermination of the series coefficients c_(n). Specifically, afterw({right arrow over (x)}) is chosen, c_(n) can be determined by solvinga set of linear equations in the following form:

Hc=d  (2)

[0025] where c and d are N-element vectors containing respectively theseries of coefficients c_(n) and the measured data d_(n) with N beingthe total number of measured data points; and H is an N×N matrix (with aToeplitz or block Toeplitz structure) constructed from the Fouriertransform samples of the weighting function w({right arrow over (x)}).For k-space data truncated only along one direction, Eq. (2) can besolved efficiently because H is relatively small (say 64×64) and has aHermitian Toeplitz structure. For data sets truncated in multipletime-dependent directions (e.g. in 3D imaging), H is significantlylarger (say, 4096×4096) and the Toeplitz structure is also lost; as aresult, solving Eq. (2) directly could present a significantcomputational problem for online processing using computers in an MRIscanner. Building on the work in Madore and Pelc, Proceedings of 8thAnnual Meeting of ISMRM, Denver, 2000, p. 297, we will first describe afast algorithm for finding an approximate solution to Eq. (2) and thenuse it for GS model-based image reconstruction.

[0026] First, we rewrite the GS model in the following form:

I({right arrow over (x)})=w({right arrow over (x)})C({right arrow over(x)})  (3)

[0027] Where${C\left( \overset{\rightarrow}{x} \right)} = {\sum\limits_{n}{c_{n}^{{- }\quad 2\quad \pi \quad {{\overset{\rightarrow}{k}}_{n} \cdot \overset{\rightarrow}{x}}}}}$

[0028] is the contrast modulation function, which is to be updated fromone frame to another. Let Î({right arrow over (x)})=I({right arrow over(x)})*h({right arrow over (x)}) be the low-resolution Fourierreconstruction of I({right arrow over (x)}) from the measured k-spacedata, where h({right arrow over (x)}) is the point spread function ofthe data truncation window. We have

{circumflex over (I)}({right arrow over (x)})=[w({right arrow over(x)})C({right arrow over (x)})]*h({right arrow over (x)})  (4)

[0029] Invoking the inherent assumption in the GS model that C({rightarrow over (x)}) is a smooth function [in the sense that C({right arrowover (x)}) does not vary significantly over any region of the shape andsize of the effective spatial support of h({right arrow over (x)})], wecan factor C({right arrow over (x)}) out of the convolution integral inEq. (4). Consequently, we have

{circumflex over (I)}({right arrow over (x)})≈[w({right arrow over(x)})*h({right arrow over (x)})]C({right arrow over (x)})={circumflexover (w)}({right arrow over (x)})C({right arrow over (x)})  (5)

[0030] where w({right arrow over (x)})=w({right arrow over(x)})*h({right arrow over (x)}) is a low-resolution version of w({rightarrow over (x)}), which is obtained by discarding the high spatialfrequency data that produced w({right arrow over (x)}). From Eq. (5), weimmediately obtain $\begin{matrix}{{C\left( \overset{\rightarrow}{x} \right)} = \frac{\hat{I}\left( \overset{\rightarrow}{x} \right)}{\hat{w}\left( \overset{\rightarrow}{x} \right)}} & (6)\end{matrix}$

[0031] The above equation, first proposed by Madore and Pelc, supra,provides a fast way to evaluate C({right arrow over (x)}) because themain computation here is the two FFTs used to calculate Î({right arrowover (x)}) and ŵ({right arrow over (x)}). In fact, for processing serialdata sets, ŵ({right arrow over (x)}) needs to be computed only oncealthough Î({right arrow over (x)}) has to be updated for each dataframe.

[0032] The solution provided by Eq. (6) could be unstable due todivision by zero (or numbers close to be zero). In accordance with afeature of the invention a straightforward way to deal with the problemis to add a “regularization” constant λ to the magnitude of ŵ({rightarrow over (x)}) such that $\begin{matrix}{{C\left( \overset{\rightarrow}{x} \right)} = \frac{\hat{I}\left( \overset{\rightarrow}{x} \right)}{\left\lbrack {{{\hat{w}\left( \overset{\rightarrow}{x} \right)}} + \lambda} \right\rbrack ^{\quad \theta \quad {(\overset{\rightarrow}{x})}}}} & (7)\end{matrix}$

[0033] where θ({right arrow over (x)}) is the phase function of ŵ({rightarrow over (x)}). In practice, only the magnitude component of C({rightarrow over (x)}) is used (see Step (a) of the fast RIGR algorithm of theinvention, as described below); we may, therefore, directly evaluateC({right arrow over (x)}) as $\begin{matrix}{{C\left( \overset{\rightarrow}{x} \right)} = \frac{{\hat{I}\left( \overset{\rightarrow}{x} \right)}}{{{\hat{w}\left( \overset{\rightarrow}{x} \right)}} + \lambda}} & (8)\end{matrix}$

[0034] Clearly, a large λ will give a more stable C({right arrow over(x)}), but at the expense of accuracy. We use a small λ (say, λ=10⁻²μwith μ being the mean value of |ŵ({right arrow over (x)})|. Ifnecessary, a median filter can be applied afterward to remove any noisespikes resulting from small values in the denominator of Eq. (8). Someadditional features of the invention have proved very useful for furtherimproving the estimate of C({right arrow over (x)}) (and avoidingspurious features in the resulting reconstruction). First, as in theoriginal RIGR algorithm (see Liang and Lauterbur, IEEE Trans Med Imaging1994; 13:677-686), it is desirable to choose a nonnegative function forw({right arrow over (x)}). Second, to reduce Gibbs' ringing in Î({rightarrow over (x)}) and ŵ({right arrow over (x)}), the Fourier data shouldbe weighted with a window function (e.g., a mild Hamming windowH(k)=0.8+0.2 cos (πk/k_(max)) in the ID case, where it is assumed thatmeasured data are available for −k_(max)≦k≦k_(max)) before the inverseFFT is applied.

[0035] After C({right arrow over (x)}) is determined from Eq. (7), a GSmodel-based reconstruction can be generated from Eq. (3) directly.However, the reconstruction as calculated will not be data-consistentbecause of the approximation, regularization, and filtering used in theprocess. It is therefore often desirable to reinforce thedata-consistency constraints by using the model-based reconstructiononly for data extrapolation; that is, in this embodiment we use themodel function in Eq. (3) to generate the unmeasured high-frequency datafor a given dynamic data set so that a “full” data set is available forFourier processing. In addition to assuring data consistency for theresulting reconstruction, the data extrapolation step will also allow usto reintroduce any valid phase constraints from the reference data tofurther improve resolution, another important feature of the invention.Based on the above discussion, a basic algorithm for GS model-basedreconstruction and two extensions are presented below, in accordancewith the invention.

[0036] A fast RIGR algorithm in accordance with the invention isapplicable to the situation where a single high-resolution imageI_(r)({right arrow over (x)}) and a sequence of reduced k-space datasets: d_(q)({right arrow over (k)}) for 1≦q≦Q, are collected. Thealgorithm requires estimation of the GS model function, dataextrapolation, and FFT reconstruction, as follows:

[0037] (a) Estimation of the GS model function:

[0038] Set w({right arrow over (x)})=|I_(r)({right arrow over (x)})|;

[0039] Calculate C({right arrow over (x)}) according to Eq. (8).

[0040] Form a model-based reconstruction I_(gs)({right arrow over (x)})according to Eq. (3).

[0041] (Optional) Apply a median filter (with a 3×3 mask) to C({rightarrow over (x)}) or I_(gs)({right arrow over (x)}).

[0042] (b) Data extrapolation:

[0043] Introduce the reference phase θ_(r)({right arrow over (x)}) toI_(gs)({right arrow over (x)}) such that Î_(q)({right arrow over(x)})=I_(gs)({right arrow over (x)})e^(iθ) _(r)^(({right arrow over (x)})).

[0044] Generate a full k-space data set {circumflex over (d)}_(q)({rightarrow over (k)}) from Î_(q)({right arrow over (x)}).

[0045] Replace {circumflex over (d)}_(q)({right arrow over (k)}) withd_(q)({right arrow over (k)}) at the k-space locations whered_(q)({right arrow over (k)}) is measured.

[0046] (c) FFT reconstruction:

[0047] Apply the inverse FFT to the merged data set {circumflex over(d)}_(q)({right arrow over (k)}) obtained in Step (b) to generate thefinal reconstruction.

[0048] In some dynamic imaging applications, it is possible to acquiretwo (or more) high-resolution images at different times of the dynamicimaging process. For example, in contrast-enhanced MRI, one can acquirea high-resolution pre-contrast image as the baseline reference andanother high-resolution reference image after the rapid dynamic wash-inphase is completed. To reconstruct the dynamic images between the tworeferences, two algorithms, known as TRIGR (Two-reference RIGR, (seeHanson et al., Magn Reson Med 1996; 36:172-175) and TRIGR+ (see Liang,Proceedings of Workshop on Minimum MR Data Acquisition Methods, MarcoIsland, Fla. 2001; 23-26), have been proposed. Based on the fast RIGRalgorithm of the invention, we can reformulate the TRIGR and TRIGR+algorithms as follows.

[0049] A fast TRIGR algorithm is applicable to the situation where twohigh-resolution references I_(r,1)({right arrow over (x)}) andI_(r,2)({right arrow over (x)}), and a sequence of reduced k-space datasets: d_(q)({right arrow over (k)}) for 1≦q≦Q, are collected. It isfurther assumed that I_(r,1)({right arrow over (x)}) can be treated asthe baseline such that I_(r,2)({right arrow over (x)})−I_(r,1)({rightarrow over (x)}) highlights areas of dynamic signal variations. Thisassumption is often valid for contrast-enhanced MRI applications [8,11].

[0050] (a) Set Ĩ_(r)({right arrow over (x)})=I_(r,2)({right arrow over(x)})−I_(r,1)({right arrow over (x)}).

[0051] (b) Remove the baseline signal from the dynamic data sets: {tildeover (d)}_(q)({right arrow over (k)})=d_(q)({right arrow over(k)})−d_(r,1)({right arrow over (k)}) for the measured {right arrow over(k)} values and 1≦q≦Q.

[0052] (c) Feed Ĩ_(r)({right arrow over (x)}) and {tilde over(d)}_(q)({right arrow over (k)}) to the fast RIGR algorithm to obtain aGS model-based reconstruction.

[0053] (d) Add I_(r,1)({right arrow over (x)}) to the image obtained in(c) to generate the final reconstruction.

[0054] The input to a fast TRIGR+ algorithm is the same as that to theTRIGR algorithm. However, here we no longer assume that one of referenceimages is a baseline scan. This algorithm is modified from the originalTRIGR+ algorithm (see L. Cong, Supra) and can be applied to data withlarge signal variations from one image to another, such as indiffusion-weighted imaging with multiple diffusion weightings.

[0055] (a) Feed I_(r,1)({right arrow over (x)}) and d_(q)({right arrowover (k)}) to the fast RIGR algorithm to obtain I_(q,1)({right arrowover (x)}).

[0056] (b) Feed I_(r,2)({right arrow over (x)}) and d_(q)({right arrowover (k)}) to the fast RIGR algorithm to obtain I_(q,2)({right arrowover (x)}).

[0057] (c) Set I_(q)({right arrow over (x)})=αI_(q,1)({right arrow over(x)})+(1−α)I_(q,2)({right arrow over (x)}) where α=e₂/(e₁ 30 e₂) withe₁=∥d_(q)({right arrow over (k)})−d_(r,1)({right arrow over (k)})∥₂ ande₂=∥d_(q)({right arrow over (k)})−d_(r,2)({right arrow over (k)})∥₂.Here it is understood that ∥·∥₂ is defined over the set of {right arrowover (k)} values at which d_(q)({right arrow over (k)}) is measured.

[0058] It is useful to state some known results before we analyze thecomputational complexity of the proposed algorithms of the invention:

[0059] The computational complexity of an N-point FFT is O(N log N).

[0060] The computational complexity of solving a Hermitian Toeplitzsystem of N simultaneous equations in N unknowns is O(N log N) [TheLevinson algorithm has a complexity of O(N²)].

[0061] The computational complexity of solving an arbitrary system of Nsimultaneous equations in N unknowns is O(N³) [or more precisely,O(N^(p)) for 2<p≦3].

[0062] The proposed fast RIGR algorithm (and its extensions) is anFFT-based algorithm, and it is easy to show that its computationalcomplexity is O(N log N) for processing a single 1D single with N beingthe number of reference encodings. For processing a 2D data set with Mcolumns, the overall computational complexity is O(MN log N) or O(N² logN) if M=N. Therefore, the computational complexity of the proposed fastRIGR algorithm is essentially the same as that of the Keyhole algorithm,although Keyhole may require fewer FFTs than the proposed fast RIGRalgorithm. The same can be said about the fast TRIGR and TRIGR+algorithms.

[0063] The original RIGR algorithm (and its variants) involves solvingEq. (2). In the case of ID data truncation, the coefficient matrix H hasa Hermitian Toeplitz structure, and Eq. (2) can be solved in O(L log L)operations or O(L²) operations if the Levinson algorithm is used [14],where L is the number of dynamic encodings. With the solution of Eq.(2), O(N log N) operations are needed to evaluate I({right arrow over(x)}) in Eq. (3), where N is the number of reference encodings. Notingthat N>L in this case, the overall computational complexity of theoriginal RIGR algorithm can be expressed simply as O(N log N), or O(MNlog N), for processing a 2D data set of M columns, which is similar tothat of the proposed fast RIGR algorithm or the Keyhole algorithm.Therefore, there is no significant computational advantage (or penalty)of the proposed fast RIGR algorithm over the original RIGR algorithm fordata sets with ID truncation.

[0064] The situation, however, is drastically different for processingdata sets truncated in multiple directions. For example, for a 2Ddynamic data set with L_(x)×L_(y) encodings, H is an L×L matrix withL=L_(x)L_(y). Although H has a block Toeplitz structure, Eq. (2) isoften solved using a general linear system solver and, therefore, thecomputational complexity of the RIGR algorithm is, in this case,approximately O(L³). In contrast, the computational complexity of theproposed algorithm(s) is O(N² log N) assuming that the reference dataset has N×N encodings. Because L is on the order of or much larger thanN, the reduction in computational complexity by the proposedalgorithm(s) is very significant; for example, for L=16×16=256 and${N = 256},{\frac{L^{3}}{N^{2}\quad \log \quad N} = 32.}$

[0065] The running time of an algorithm is dependent on itsimplementation and the computing platform used. Nevertheless, with thesignificantly reduced computational complexity, the proposed fast RIGRalgorithm (and its extensions) are expected to offer a significantspeedup in GS model-based image reconstruction over the originalalgorithms for data sets truncated in multiple directions. We haveobserved up to a factor of ten speed-ups by the proposed algorithms inrunning our unoptimized Matlab implementation of the algorithms.

[0066] Following are some representative results from practical imagingdata to demonstrate the effectiveness of the algorithms in accordancewith the invention.

[0067] The first data set was acquired using a fast gradient-echosequence following intravenous administration of a contrast agent. Westarted with a full data set consisting of 20 temporal data frames, eachwith 144 phase-encoded echo signals. As the echoes are asymmetric (with26 and 128 data points respectively before and after the echo peak), weused a projection-onto-convex-sets algorithm for partial-Fourierreconstruction along the readout direction (the vertical direction inFIG. 2). To test the proposed algorithms, a large portion of theacquired data was discarded, and two high-resolution images (top row ofFIG. 2), were used as the reference, of which the first is apre-injection image and the second is a post-injection image acquired atthe end of the scan. FIGS. 2b-e shows a set of reconstructions fromseveral algorithms at three time points during transit of the contrast.The Fourier reconstructions from the full data set (FIG. 2b) are alsoincluded as the gold standard. As can be seen from FIG. 2c, the Keyholereconstructions with 16 phase encodings contain noticeablereconstruction errors (over-enhancements of blood vessels in the firstimage and under-enhancements in the second image). These artifacts weresignificantly reduced in the GS model-based reconstructions produced byboth the proposed fast TRIGR algorithm (FIG. 2d) and the original TRIGRalgorithm (FIG. 2e). In fact, there is no noticeable difference in imagequality between FIGS. 2d and e. This observation is consistent with theresults shown in FIG. 3, which indicates that the proposed fast TRIGRalgorithm performs as well as the original TRIGR algorithm for differentlevels of data truncation of this data set.

[0068] The second data set was obtained from a diffusion-weightedimaging experiment, where “dynamic” signal changes were caused byvarying the amplitude of a diffusion-weighting gradient. This data setwas chosen to test how well the algorithms handle large signalvariations with respect to the reference(s).

[0069]FIG. 4a shows two high-resolution references, one acquired with aminimum diffusion weighting and the other with a heavy diffusionweighting. FIGS. 4c-e shows the reconstructions of threediffusion-weighted images from 16 phase-encoded signals using Keyholeand the GS model-based methods. For comparison, the correspondinghigh-resolution Fourier images reconstructed from 128 phase-encodedsignals are also included (FIG. 4b). As can be seen, the Keyholereconstructions (FIG. 4c) show significant image blurring because oflarge “dynamic” signal variations from the reference [the second imagein (3 a)] Keyhole reconstructions using the first reference are notshown here as they contain stronger artifacts. This blurring artifactwas significantly reduced in the GS model-based reconstructions producedby both the original and fast TRIGR+ algorithms (FIGS. 4d-e). FIG. 5shows the reconstruction error curves calculated from the entirediffusion-weighted data set with 6 frames. As expected, the TRIGRalgorithms have smaller reconstruction errors than Keyhole. Note alsothat the reconstruction error of the fast TRIGR+ algorithm is slightlyhigher than that of the original TRIGR+ algorithm because the smoothnessassumption made in Eq. (5) is violated in this application. The fastalgorithms of the invention would seem to be more ideal incontrast-enhanced applications (e.g. FIG. 2), where the contrastmodulation function C(x) is intrinsically smooth and the assumption madein Eq. (5) holds well, than in the diffusion application (FIG. 4).

[0070] The GS model is an optimal way (according to the minimumcross-entropy principle) to solve the image reconstruction problem indata-sharing Fourier imaging. The invention includes three related fastalgorithms for calculating GS model-based reconstructions. Experimentalresults indicate that the proposed algorithms can produce images ofhigher quality than those from the Keyhole algorithm and comparable tothose by the original GS model-based algorithms. The proposedalgorithms, having the same computational complexity as the Keyholealgorithm make the GS model more practical for processing k-space datatruncated in multiple directions.

[0071] While the invention has been described with reference to specificapplications, the description is illustrative of the invention and isnot limiting the invention. Various applications and modifications mayoccur to those skilled in the art without departing from the true scopeand spirit of the invention as defined by the appended claims.

What is claimed is:
 1. A method of image reconstruction for fast dynamicimaging using a generalized-series (GS) model comprising the steps of:a) acquiring data for at least one high resolution reference image,I_(r,q)({right arrow over (x)}), q≧1. b) acquiring one or more k-spacedata sets: d_(q)({right arrow over (k)}) for 1≦q≦Q c) determining aweighting function w({right arrow over (x)}) based on the referenceimage(s). d) determining a contrast modulation function, C({right arrowover (x)}), for each data frame to absorb dynamic signal variations, ande) multiplying the weighting function w({right arrow over (x)}) by thecontrast modulation function, C({right arrow over (x)}), to obtain aninitial estimate of the image for each time point. f) using the initialreconstruction obtained in step e) and the phase information from thereference image(s) to create a full k-space data set. g) Replacing theestimated k-space data obtained in step f) with the measured data at thek-space locations where d_(q)({right arrow over (k)}) is measured toachieve data consistency. h) Fourier transforming the merged data setobtained in step g) to generate the desired image function.
 2. Themethod as defined by claim 1 wherein in step a) the reference data arecollected using the standard phase-encoding or frequency-encodingmethods, and the number of encodings is determined by the desiredspatial resolution.
 3. The method as defined by claim 1 wherein in stepb) the dynamic data sets are collected using the standard phase-encodingor frequency-encoding methods, spiral scanning methods, and the numberof encodings is determined by the desired temporal resolution of theexperiment.
 4. The method as defined by claim 1 wherein in step c) theweighting function is taken as the magnitude of the reference image oras the magnitude of the difference between two references.
 5. The methodas defined by claim 1 wherein in step d) the contrast modulationfunction is given by${{C\left( \overset{\rightarrow}{x} \right)} = \frac{{\hat{I}\left( \overset{\rightarrow}{x} \right)}}{{\overset{\bullet}{w}\left( \overset{\rightarrow}{x} \right)}}},$

with a) Î({right arrow over (x)}) being the low-resolution Fourierreconstruction of I_(q)({right arrow over (x)}) from the measuredk-space data d_(q)({right arrow over (k)}); and b) ŵ({right arrow over(x)}) being the low-resolution version of w({right arrow over (x)}),which is obtained by discarding the high spatial frequency data thatproduced w({right arrow over (x)}).
 6. The method as defined by claim 1wherein in step d) the contrast modulation function is given by${{C\left( \overset{\rightarrow}{x} \right)} = \frac{{\hat{I}\left( \overset{\rightarrow}{x} \right)}}{{{\hat{w}\left( \overset{\rightarrow}{x} \right)}} + \lambda}},$

with a) Î({right arrow over (x)}) being the low-resolution Fourierreconstruction of I_(q)({right arrow over (x)}) from the measuredk-space data d_(q)({right arrow over (k)}); b) ŵ({right arrow over (x)})being the low-resolution version of w({right arrow over (x)}), which isobtained by discarding the high spatial frequency data that producedw({right arrow over (x)}); and c) λ being a regularization constant. 7.The method as defined by claim 6 wherein λ being a regularizationconstant whose value is roughly 10⁻²μ where μ is the mean value of|ŵ({right arrow over (x)})|.
 8. The method as defined by claim 1 whereinin step a) two high resolution images, I_(r,1)^(({right arrow over (x)})) and I_(r,2) ^(({right arrow over (x)})), areacquired at different times with I_(r,1) ^(({right arrow over (x)}))being a base-line image; wherein in step c) the weighting function ischosen to be |I_(r,2) ^(({right arrow over (x)}))−I_(r,1)^(({right arrow over (x)}))|; in step d) the contrast modulationfunction is determined from {right arrow over (d)}_(q)({right arrow over(k)})=d_(q)({right arrow over (k)})−d_(r,1)({right arrow over (k)}), andin step h) I_(r,1) ^(({right arrow over (x)})) is added to thereconstruction of the fast RIGR algorithm to produce the final image foreach temporal frame.
 9. The method as defined by claim 1 wherein in stepa) two high resolution images, I_(r,1) ^(({right arrow over (x)})) andI_(r,2) ^(({right arrow over (x)})) are acquired at different times; twoGS model-based reconstructions are produced by the fast RIGR algorithmof the invention using |I_(r,1) ^(({right arrow over (x)}))| and|I_(r,2) ^(({right arrow over (x)}))|, respectively, as the weightingfunction, and the final reconstruction is determined as a weighted sumof the two GS model-based reconstructions.
 10. The method as defined byclaim 9 where the weighting coefficients are determined according to∥d_(q)({right arrow over (k)})−d_(r,1)({right arrow over (k)})∥₂ and∥d_(q)({right arrow over (k)})−d_(r,2)({right arrow over (k)})∥₂.